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Abstract 

We present a meshless simulation method for multiphase fluid-particle flows 
coupling Smoothed Particle Hydrodynamics (SPH) and the Discrete Element 
Method (DEM). Rather than fully resolving the interstitial fluid, which is often 
infeasible, the unresolved fluid model is based on the locally averaged Navier 
Stokes equations, which are coupled with a DEM model for the solid phase. 
In contrast to similar mesh-based Discrete Particle Methods (DPMs), this is a 
purely particle-based method and enjoys the flexibility that comes from the lack 
of a prescribed mesh. It is suitable for problems such as free surface flow or flow 
around complex, moving and/or intermeshed geometries. It can be used for both 
one and two-way coupling and is applicable to both dilute and dense particle 
flows. A comprehensive validation procedure for fluid-particle simulations is 
presented and applied to the SPH-DEM method, using simulations of single 
and multiple particle sedimentation in a 3D fluid column and comparison with 
analytical models. Millimetre-sized particles are used along with three different 
test fluids: air, water and a water-glycerol solution. The velocity evolution for 
the single particle compares well (less than 2% error) with the analytical solution 
as long as the fluid resolution is coarser than 2 times the particle diameter. 
The multiple particle sedimentation problems (sedimentation of a homogeneous 
porous block and a Rayleigh Taylor instability) also reproduce the expected 
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terminal velocity well for porosities 0.5 < e < 1.0, but although care should be 
taken in the presence of high porosity gradients. Overall the SPH-DEM method 
successfully reproduces the expected behaviour in the sedimentation test cases, 
and promises to be a flexible and accurate tool for other fluid-particle system 
simulations. 

Keywords: SPH, DEM, Fluid-particle flow, Discrete Particle Model, 
Sedimentation, Rayleigh- Taylor instability 



1. Introduction 



Fluid-particle systems are ubiquitous in nature and industry. Sediment 
transport and erosion are important in many environmental studies and the in- 
teraction between particles and interstitial fluid affects the rheology of avalanches, 
slurry flows and soils. In industry, the efficiency of a fluidised bed process (e.g. 
Fluidized Catalytic Cracking) is completely determined by the complex two- 
way interaction between the injected gas flow and the solid granular material. 
Also, the dispersion of solid particles in a fluid is of broad industrial relevance 
to the food, chemical and painting industries, which involves in most cases three 
phases: a granular medium, the air initially present in its pores and an injected 
liquid. 

The length-scale of interest determines the method of simulation for fluid- 
particle systems. For very small scale proces ses it is feasible t o fully resolve 
the in t erstitial fluid be t ween the particles (see 



( 20101) 



Potapovl (|2001l) 



Wachmann et al 



Zhu et al 



(1999); 



Pereira et al 



(|1998l ) for a few examples of particle 



or pore-scale simulations). However, for many applications the dynamics of in- 
terest occur over length scales much larger than the particle diameter and the 
computational effort required to resolve the pore-scale is too great. It then 
becomes necessary to use unresolved, or mesoscale, fluid simulations. This 
mesoscale is the focus of this paper and the domain of applicability for the 
SPH-DEM method. At even larger length scales of interest (macroscale) it be- 
comes infeasible to model the granular material as a discrete collection of grains 
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and instead a continuum model is used in a two-fluid model. However, it must 
be noted that while this approach might be computationally necessary in many 
cases, it can fail for some systems involving dense granular flow, where existing 
continuum models for granular material do not adequately reproduce impor- 
tant material properties such as anisotropy, history dependency, jamming and 
segregation. 

Fluid-particle simulations at the mesoscale are often given the term Discrete 
Particle Models (DPM). These models fully resolve the individual solid parti- 
cles using a Lagrangian model for the solid phase. The fluid phase does not 
resolve the interstitial fluid, but instead models the locally averaged Navier- 
Stokes equations and is coupled to the solid particles using appropriate drag 
closures. Most of the prior work on DPMs have been done using grid-based 



pers by 


Tsuii et al. 


(1993 


), 2 


(2000 


) or 


Chu and Yu 


(2008) 



Xu (1997 



2000), 



Hoomansl (1996); 



Hoomans et al 



Fixed pore flow simulations (where the geometry of the solid particles is 
unchangin g over time) us ing SP H for the (unreso lved) fluid phase have been de- 



scribed by 



Li et al 



(|2007j ) and l Jiang et al.1 (120071). but the se do not allow for the 



motio n and collision of solid grains. 



Clearv et al 



hOO®) and 



Fernandez et al 



( 20111 ) simulate slurry flow at the mesoscale using SPH and DEM in SAG mills 
and through industrial banana screens, but only perform a one-way coupling 
between the solid and fluid phases. 

The DPM model presented in this paper is based on the locally averaged 
Navier-Stokes (A VNS) equations that were first derived by Anderson and Jack- 
son in the sixties ([Anderson and JacksonL 119671 ) , and have been used with great 
success to mod el the comp l ex flu id-particle interactions occurring in industrial 



fluidizcd beds (Deenetal 



20071) . Anderson and Jackson defined a smoothing 



operator identical to that used in SPH and used it to reformulate the NS equa- 
tions in terms of smoothed variables and a local porosity field (porosity refers 
to the fraction of fluid in a given volume) . Given its theoretical basis in kernel 
interpolation, it is natural to consider the use of the SPH method to solve the 
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AVNS equations, coupled with a DEM model for the solid phase. 

The coupling of SPH and DEM results in a purely particle-based solution 
method and therefore enjoys the flexibility that is inherent in these methods. 
This is the primary advantage of this method over existing grid-based DPMs. 
In particular, the model described in this paper is well suited for applications 
involving a free surface, including (but not limited to) debris flows, avalanches, 
landslides, sediment transport or erosion in rivers and beaches, slurry transport 
in industrial processes (e.g. SAG mills) and liquid-powder dispersion and mixing 
in the food processing industry. 

Another advantage of using a DPM, or mesoscale simulation, is of course 
the reduced computational requirements over a fully resolved simulation. For 
example, we have found that in general a fluid resolution of h = 2d minimises 
the error in the SPH-DEM method, where d is the solid particle diameter. For 
a fully resolved simulation the interstitial fluid must be resolved, and therefore 
the fluid resolution would need to be at least h = 0.2d, which scales the number 
of computational nodes (for the fluid) by a factor of 1000. 

Figure [1] shows a SPH-DEM simulation applied to a liquid-powder mixing 
problem in the food processing industry, taken from a simulation of a water jet 
injected in a granular bed whose pores are initially filled with air. To predict 
the shape of the front correctly, one has to consider the free surface and the 
absence of dissipation on the air side, both included in the SPH-DEM model. 
Even more complex (realistic) injection geometries are easily incorporated into 
the simulation with no additional effort. Moreover, using DEM enables studying 
the effect (on the initial liquid front propagation) of packing and top surface 
inhomogeneities that can be generated during pouring, unlike simpler "porous 
media" -like approaches. Polydispersity can also be included by altering the 
radius of the simulated g rains and using an applicable drag term (e.g. see 
Van der Hoef et al.1 1200§)) 

Sections describes the AVNS equations and the SPH and DEM models 
for the fluid and solid phases and the coupling between them. Section [4] in- 
troduces the test cases, Section [5] describes the results for the Single Particle 
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Figure 1: Example of a two-phase SPH-DEM simulation of a water jet (bottom) injected into 
a granular bed. 



Sedimentation test case, Section [6] the results for Multiple Particle Sedimen- 
tation and Section [7] describes the simulation of a Rayleigh Taylor Instability 
using solid particles sedimenting into a clear fluid. 

2. Governing Equations 

2.1. The Locally Averaged Navier-Stokes Equations 

Here we describe the governing equati ons for the fluid phase, the loc ally aver- 



aged Navier-Stokes equations derived bv I Anderson and Jackson! (|1967l K Ander 



son and Jackson defined a local averaging based on a radial smoothing function 
g(r). The function g(r) is greater than zero for all r and decreases monotonically 
with increasing r, it possesses derivatives g n (r) of all orders and is normalised 
so that / g(r)dV = 1. 

The local average of any field a' defined over the fluid domain can be obtained 
by convolution with the smoothing function 
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e(x)a(x) = / a'(y)g(x - y)dV y , (1) 
J Vf 

where x and y are position coordinates (here one dimensional for simplicity). 
The integral is taken over the volume of interstitial fluid Vf and e(x) is the 
porosity. 

e(x) = l- f g(x~y)dV v , (2) 
Jv s 

where V s is the volume of the solid particles. 

In a similar fashion, the local average of any field a'(x) defined over the solid 
domain is given by 



(l-e(x))a(x) = / a'(y)g(x-y)dV y , (3) 
Jv s 

where the integral is taken over the volume of the solid particle^ 



Applying this averaging method to the Navier-Stokes equations. 



Anderson and Jackson 



(Il967t ) derived the following continuity equation in terms of locally averaged 



variables 

d(ep f ) 



V • (ep f u) = 0, (4) 



Of 

where pf is the fluid mass density and u is the fluid velocity. 
The corresponding momentum equation is 

ep f \j£ + u • Vuj = -VP + V • r - nf + ep/g, (5) 

where P is the fluid pressure, r is the viscous stress tensor and nf is the 
fluid-particle coupling term. We use a Newtonian fluid where r = /iV • u. We 
neglect Reynolds-like terms and do not consider turbulent flow. The coefficient 
for the coupling term n is the local average of the number of particles per unit 
volume and f is the local mean value of the force exerted on the particles by the 
fluid. This force includes all effects, both static and dynamic, of the particles 
on the fluid, the details of which can be seen in Eq. ([23]) . 
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2.2. Smoothed Particle Hydrodynamics 



Smoothed Par ticle Hydrodynamics (jGingold and Monaghan , 



Monaghan , 



1977 



Lucy 



1977: 



2005) is a Lagrangian scheme, whereby the fluid is discretised into 
"particles" that move with the local fluid velocity. Each particle is assigned a 
mass and can be thought of as the same volume of fluid over time. The fluid 
variables and the equations of fluid dynamics are interpolated over each par- 
ticle and its nearest neighbours using a smoothing kernel W{r, h), where h is 
the smoothing length scale. Like g(r) in the AVNS equations, the SPH ker- 
nel is a radial function that decreases monotonically and is normalised so that 
/ W{r,h)dV = 1. 

Unlike g{r) and to reduce the computational burden of the method, the 
SPH kernel is normally defined with a compact support and a finite number of 
derivatives. 

In SPH, a fluid variable ^4(r) (such as momentum or density) is interpolated 
using the kernel W 

A(r) = J A(r')W(r-r',h)dr'. (6) 

To apply this to the discrete SPH particles, the integral is replaced by a sum 
over all particles, commonly known as the summation interpolant. To estimate 
the value of the function A at the location of particle a (denoted as A a ), the 
summation interpolant becomes 

A a = £m 6 — W„ 6 (ft«)> (7) 
b Pb 

where nib and pb are the mass and density of particle b. The volume element 
dr' of Eq. (|6]) has been replaced by the volume of particle b (approximated by 
^■), equivalent to the normal trapezoidal quadrature rule. The kernel function 
is denoted by W a b(h) = W(r a — rf,,/i). The dependence of the kernel on the 
difference in particle positions is not explicitly stated for readability. Due to 
the limited support of W, particle neighbourhood search methods as standard 
in SPH or DEM can be applied to optimize the summation in Eq. ([6]). 
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The accuracy of the SPH interpolant depends on the particle positions within 
the radius of the kernel. If there is not a homogeneous distribution of particles 
around particle a (for example, it is on a free surface), then the interpolation 
can be compromised. 



The interpolation can be improved by using a Shepard correction ([Shepardl . 



1968J), originally devised as a low cost improvement to data fitting. This cor- 
rection divides the interpolant by the sum of kernel values at the SPH particle 
positions, so the summation interpolant becomes 

This correction ensures that a constant field will always be interpolated 
exactly, and improves the interpolation accuracy of other, non-constant fields. 

3. SPH-DEM Model 

3.1. SPH implementation of the AVNS equations 

SPH is based on a similar local averaging technique as the AVNS equations, 
so it is natural to convert the interpolation integrals in Eqs. (TT|) and ^ to SPH 
sums using a smoothing kernel W(r, h) in place of g(r). 

To calculate the porosity e a at the center position of SPH/DEM particle a, 
the integral in Eq. ^ is converted into a sum over all DEM particles within 
the kernel radius and becomes 



e a = l-J2w a] (h c )V ] , (9) 

3 

where Vj is the volume of DEM particle j. For readability, sums over SPH 
particles use the subscript b, while sums over surrounding DEM particles use 
the subscript j. Note that we have used a coupling smoothing length h c to 
evaluate the porosity, which sets the length scale for the coupling terms between 
the phases. Here we set h c to be equal to the SPH smoothing length, but in 
practice this can be set within a range such that h c is large enough that the 
porosity field is smooth but small enough to resolve the important features of 
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the porosity field. For more details on this point please consult the numerical 
results of the test cases and the conclusions of this paper. 

Applying the local averaging method to the Navier-Stokes equations, Ander- 
son and Jackson derived the continuity and momentum equations shown in Eqs. 
([4]) and ([5]) respectively. To convert these to SPH equations, we first define a 
superficial fluid density p equal to the intrinsic fluid density scaled by the local 
porosity p = epf. 

Substituting the superficial fluid density into the averaged continuity and 
momentum equations reduces them to the normal Navier-Stokes equations. 
Therefore, our approach is to use th e stan dard weakly compressible SPH equa- 



tions, see ( Robinson and Monaghan , 



20111 ). using the superficial density for the 



SPH particle density and adding terms to model the fluid-particle drag. 

The rate of change of superfic ial density i s calculated using the variable 



smoothing length terms derived by 



Pried \20m 



(10) 



Dt n a 

b 

where the capitals in the time derivative denote the material derivative. 
Uab = u a — u b and Q, a is a correction factor due to the gradient of the smoothing 
length 



O a = 1 



dh a ^ dW ab (h a ) 



dp a 



dh a 



(11) 



Neglecting gravity, the SPH acceleration equation becomes 



du a \ - 



Pa 



VaPl 



n o6 v a w ab {h a ) + 



n Qh v a w ab {h b ) 



+f a /m a , 
(12) 

where f a is the coupling force on the SPH particle a due to the DEM particles 
(see Section 1575)1 . The viscous term H ab models the divergence of the viscous 



stress tensor in Eq. ([5]) is calculated using the term proposed by 



Monaghan 
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(|1997f) . which is based on the dissipative term in shock solutions based on Ric- 
mann solvers. For this viscosity 



n ab = -a^^- (13) 

2 Pab\ r ab\ 

where u S i g = c s + u n /\r a b\ is a signal velocity that represents the speed 
at which information propagates between the particles. The normal velocity 
difference between the two particles is given by u n = u Q b • r a f>. The constant a 
can be related to the dynamic viscosity of the fluid /i using 



jj, = pahc s /S, 



(14) 



where S — 112/15 for two dimensions and S = 10 for three (jMonaghan . 



20051) . For some of the reference fluids we have chosen to simulate in this paper 
it was found that the physical viscosity was not sufficient to stabilise the results 
(see Section 16. 3p , and it was necessary to add an artificial viscosity term with 
otart =0.1. However, this viscosity term is only applied when the SPH particles 
are approaching each other (i.e. u a b ■ r a b < 0) so that the dissipation due to the 
artificial viscosity is reduced while still stabilising the results. 

The fluid pressure in Eq. (TT2"j) is calculated using the weakly compressible 
equation of state. This equation of state defines a reference density p at which 
the pressure vanishes, which must be scaled by the local porosity to ensure that 
the pressure is constant over varying porosity. 



Pa = B 



Pa 

e a Po 



1 



(15) 



The scaling factor B, is free a-priori and is set so that the density variation 
from the local reference density is less than 1 percent, ensuring that the fluid is 
close to incompressible. For this, in terms of B, the local sound speed is 



dP 



-yB 



(16) 



P=e a pa 

and the fluct uations in density can be related to the sound speed and velocity 



of the particles (jMonaghan , 



2005) 
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M = K- (17) 

P c% 

Therefore, in order to keep these fluctuations less than 1% in a flow where 
the maximum velocity is u m and the maximum porosity is as always e m = 1, B 
is set to 

B = m *><. (18) 

7 

As the superficial density will vary according to the local porosity, care must 
be taken to update the smoothing length for all particles in order to maintain 
a sufficient number of neighbour particles. This is referred to as " variable- h" in 
this study. The smoothing length h a is calculated using 

K^) 1 ", (19, 

where d is the number of dimensions and a determines the resolution of the 
summation interpolant. The value used in all the simulation results presented 
here is a — 1.5. 

Recall that the SPH density is given by p — epf. Assuming a constant 
intrinsic fluid density pf, the smoothing length h is thus proportional to the 
local porosity h cx (l/e) 1 /^. 

3.2. Discrete Element Model (DEM) 

In DEM (also known as Molecular Dynamics), Newton's equations of mo- 
tion are integrated for each individual solid particle. Interactions between the 
particles involve explicit force expressions that are used whenever two particles 
come into contact. 

Given a DEM particle i with position r*, the equation of motion is 

m = ^2 c v + fj + m * g ' ( 20 ) 

3 

where m t is the mass of particle i, Cjj is the contact force between particles 
i and j (acting from j to i) and fj is the fluid-particle coupling force on particle 
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i. For the simulations presented below, we have used the linear spring dashpot 
contact model 



a, =-(fc*-/9£)ny, (21) 

where S is the overlap between the two particles (positive when the particles 
are overlapping, zero when they are not) and is the unit normal vector 
pointing from j to i.. The simulation timestep is calculated based on a typical 



contact duration t c and is given by At = with t c — n/ yj (2k/rrii) — (3/rrii. 
The timestep for the SPH method is set by a CFL condition 

5h < mm ( 0.6— ) , (22) 
a \ u sig J 

where the minimum is taken over all the particles. This is normally much 
larger than the DEM contact time, so the DEM timestep usually sets the min- 
imum timestep for the SPH-DEM method. 

See Table [T] for all the parameters and time-scales used in these simulations. 

3.3. Fluid- Particle Coupling Forces 



The force on each solid particle by the fluid is ([Anderson and Jackson , 



1967) 



f i = V i (-VP + V-T) i +f d {e i ,u s ), (23) 

where Vi is the volume of particle i. The first two terms models the effect of 
the resolved fluid forces (buoyancy and shear-stress) on the particle. For a fluid 
in hydrostatic equilibrium, the pressure gradient will reduce to the buoyancy 
force on the particle. The divergence of the shear stress is included for com- 
pleteness and ensures that the movement of a neutrally buoyant particle will 
follow the fluid streamlines. For the simulations considered in this paper this 
term will not be significant. 

The force is a particle drag force that depends on the local porosity e» 
and the superficial velocity u s (defined in the following section). This force 
models the drag effects of the unresolved fluctuations in the fluid variables and 
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is normally denned using both theoretical arguments and fits to experimental 
data. For a single particle in 3D creeping flow this term would be the standard 
Stokes drag force. For higher Reynolds numbers and multiple particle inter- 
actions this term is deter mined using fits to numerical or experimental data 



( Van der Hoef et al 



20051 ) . See Section [3^1 for further details. 
The pressure gradient and the divergence of the stress ten sor are evaluated at 
each solid particle using a Shepard corrected (|Shepardl ll968) SPH interpolation. 
Using the already given SPH acceleration equation, Eq. (fT2")l . this becomes 



(-VP + V-r) t = * ]Tm b fW4M, (24) 



b 



Ebjtw ab (h b ) , 

1 " • Uab) VaW ab (h a ) + ( ^ + U ab ) V a W ab {h b ) 



(25) 

In order to satisfy Newtons third law (i.e. the action = reaction principle), 
the fluid-particle coupling force on the fluid must be equal and opposite to the 
force on the solid particles. Each DEM particle is contained within multiple 
SPH interaction radii, so care must be taken to ensure that the two coupling 
forces are balanced. 

The coupling force on SPH particle a is determined by a weighted average 
of the fluid-particle coupling force on the surrounding DEM particles. The 
contribution of each DEM particle to this average is scaled by the value of the 
SPH kernel. 

fa = -— EtUw^CM, (26) 

Pa >Jj 
3 

where fj is the coupling force calculated for each DEM particle using Eq. 
(|2"3"|) . The scaling factor Sj is added to ensure that the force on the fluid phase 
exactly balances the force on the solid particles. It is given by 



^=E— w ^hc), (27) 
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where the sum is taken over all the SPH particles surrounding DEM particle 
j. For a DEM particle immersed in the fluid this will be close to unity. 

3.4- Fluid- Particle Drag Laws 

The drag force depends on the superficial velocity u s , which is proportional 
to the relative velocity between the phases. If Uf and are the fluid and particle 
velocity respectively, then the superficial velocity is defined as u s = £j(u/ — Uj). 
This term is used as the dependent variable in many drag laws as it is easily 
measured from experiment by dividing the fluid flow rate by the cross-sectional 
area. 

In the SPH-DEM model, the fluid velocity Uf used to calculate the superficial 
velocity, is found at each DEM particle position using a Shepard corrected SPH 
interpolation. The value of the porosity field at each DEM particle position 
is found in an identical way. 

The simplest drag law is the Stokes drag force 



f d = 3irfidu S7 (28) 

where d is the particle diameter. This is valid for a single particle in creeping 
flow. 



Coulson and Richardson 



([19931 ) proposed a drag law valid for a single parti- 



cle falling under the full range of particle Reynolds Numbers Re p — pf\u s \d/ fx. 

f d = ^d 2 p f \u s \ (l.84i? ep 031 + 0.293i?e° 06 ) 3 - 45 (29) 

For higher Reynolds numbers and multiple particles, the drag law can be 
generalised to 

f d = lc d f{e i )Kd 2 Pf \vL s \u s , (30) 

o 

where C d is a drag coefficient that varies with the particle Reynolds number 
Re p — pf\u s \d/ fj, 7 and /(e^) is the voidage function that models the interactions 
between multiple particles in the fluid. 
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A popular definition for the drag coefficient was proposed by iDallavalle 
( 1948 ) 



C d = 



0.63 



4.8 

./Rfk. 



(31) 



Di Felice proposed a void age funct i on ba sed on experimental data of fluid 



flow through packed spheres (|Di Felicd . 



1994) 



£ = 3.7- 0.65 exp 



(1.5 - log 10 Rep) 2 



(32) 
(33) 



Both the Stokes drag term (as the simplest reference case) and the com- 
bination of Dallavalle and Di Felice's drag terms are used in the simulations 
presented in this paper. A nother c omm only used drag term is giv en by a com- 
bination of drag terms by lErgunl (|l952i ) and IWen and Yul (| 19661 ) . For a — > 1 
this term and Di Felice's are identical (over all Re). As the porosity decreases 
both drag terms generally follow the same trend, although the Ergun and Wen 
& Yu model gives a larger drag force for dense systems. 

4. Validation Test Cases 

Three different sedimentation test cases were used to verify that SPH-DEM 
correctly models the dynamics of the two phases (fluid and solid particles) and 
their interactions. 

1. Single Particle Sedimentation (SPS) 

2. Sedimentation of a constant porosity block (CPB) 

3. Rayleigh Taylor Instability (RTI) 

These test cases were designed to test the particle-fluid coupling mechanics 
in order of increasing complexity. The first test case simply requires the correct 
calculation and integration of the drag force on the single particle, the single 
particle being too small to noticeably alter the surrounding fluid velocity. The 
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I 



Figure 2: Setup for test case SPS, single particle sedimentation in a fluid column. (Left) 
Perspective view, showing the fluid domain, the no-slip bottom boundary and the single 
spherical DEM particle. (Right) Top view, the grey area is the bottom no-slip boundary 



second requires that the drag on both phases and the displacement of fluid 
by the particles be correctly modelled for a simple velocity field and constant 
porosity. The third test case does the same but with a more complicated and 
time- varying velocity and porosity field. 

The first test case (SPS) models a single particle sedimenting in a fluid 
column under gravity. Figure [5] shows a diagram of the simulation domain. 
The water column has a height of h = 0.006m and the bottom boundary is 
constructed using Lennard- Jones repulsive particles (these particles are identical 
to those used by iMonaghan et al.l (J2003)). The boundaries in the x and y 
directions are periodic with a width of w = 0.004 m and gravity acts in the 
negative z direction. The single DEM particle is initialised at z = 0.8h. It has 
a diameter equal to d = 1 x 10~ 4 m and has a density p p = 2500 kg/m 3 . 

For the initial conditions of the simulation, the position of the DEM particle 
is fixed and the fluid is allowed to reach hydrostatic equilibrium. The particle 
is then released at t = s. 

Most fluid-particle systems of interest will involve large numbers of particles, 
and therefore the second test case (CPB) involves the sedimentation of multiple 
particles through a water column. In this layer of sedimenting particles 

is placed above a clear fluid region. Figure [3] shows the setup geometry. The 
fluid column is identical to the previous test case, but now the upper half of the 
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Figure 3: Setup for test cases CPB and RTI, multiple particle sedimentation in a fluid column. 

column is occupied by regularly distributed DEM particles on a square cubic 
lattice, with a given porosity e. The separation between adjacent DEM particles 
on the lattice is given by Ar = (V/ (1 — e)) 1 ^ 3 , where V is the (constant) particle 
volume. The diameter and density of the particles are identical to the single 
particle case. In order to maintain a constant porosity as the layer of particles 
falls, the DEM particles are restricted from moving relative to each other and 
the layer of particles falls as a block (only translation, no rotation of the layer) . 

The third test cases (RTI) uses the same simulation domain and initial con- 
ditions as CPB, but now the particles are allowed to move freely. This setup 
is similar in nature to the classical Rayleigh- Taylor (RT) instability, where a 
dense fluid is accelerated (normally via gravity) into a less dense fluid. The 
combination of particles and fluid can be modeled as a two-fluid system with 
the upper "fluid" having an effective density pd, and an effective viscosity /id, 
both higher than the properties of the fluid without particles. From this an 
expected growth rate can be calculated for the instability and compared with 
the simulated growth rate. See Section [7] for more details. 

For all three test cases, three different model fluids are used to evaluate the 
SPH-DEM model at different fluid viscosities and particle Reynolds numbers. 
The density and viscosities of these fluids correspond to the physical properties 
of air, water and a 10% glycerol-water solution. 
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4-1- Simulation Parameters and Timescales 

Table [T] shows the parameters used in the three test cases. Each column 
corresponds to a different model fluid. Where a value appears only in one 
column, this indicates that the parameter is constant for all the fluids. The 
particle Reynolds number is calculated using the expected terminal velocity of 
the single particle or multiple particle block. 

The standard Stokes law, Eq. (f28|). can be used to calculate the vertical 
speed of a single particle falling in a quiescent fluid. 

v (t) = [EE ~^ Vg (l - e -W / m ) , with constant b = Zirpd. (34) 

Since we are interested in a range of particle Reynolds numbers, not just at 
the Stokes limit, we also consider the Di Felice drag force, Eq. (|32l) . which is 
valid for higher Reynolds numbers and varying porosity (i.e. it considers the 
interaction of multiple particles). When the buoyancy and gravity force on the 
falling particle balance out the drag force, the particle is falling at its terminal 
velocity Equating these terms leads to a polynomial equation in terms of the 
particle Reynolds number at terminal velocity 

0.392i?ep + 6.048i?e p 5 + 23.04.Rep - ^Are 1+i = 0, (35) 

where £ is given in Eq. (|32|) and Ar — d 3 p(p p — p)g/ p 2 is the Archimedes 
number. The Archimedes number gives the ratio of gravitational forces to vis- 
cous forces. A high Ar means that the system is dominated by convective 
flows generated by density differences between the fluid and solid particles (e.g. 
Buoyancy, Rayleigh Taylor instabilities) . A low Ar means that viscous forces 
dominate and the system is governed by external forces only. 

Solving for Re p , one can find the expected terminal velocity using Re p — 
p\u t \d/fi. 

Note that a range of porosities is used for test cases CPB and RTI, and this 
results in a range of particle Reynolds numbers as the terminal velocity depends 
on the porosity. 



18 



Table 1: Relevant parameters and timescales for the simulations using different fluids. Parameters appearing only in one column are kept constant 
for all fluids. 





Notation 


Units 


Air 


Water 


Water + 10% Glycerol 


Box Width 


w 


m 


4 x 10" 3 






Box Height 


h 


m 


6 x 10 -3 






Fluid Density 


P 


kg/m 3 


1.1839 


1000 


1150 


Fluid Viscosity 


P- 


Pa-s 


1.86 x icr 5 


8.9 x 10" 4 


8.9 x 10" 3 


Particle Density 


Pp 


kg/m 3 


2500 






Particle Diameter 


d 


m 


1.0 x 10~ 4 






Spring Stiffness 


k 


kg/s 2 


1.0 x 10~ 4 






Spring Damping 


P 


kg/s 









Porosity 


e 




0.6-1.0 






Calculated Terminal Velocity (Eq. I35|l 


\ut\ 


m/s 


0.102-0.5 


1.3 x 10~ 3 -7.6 x 10~ 3 


1.3 x 10~ 4 -8.4 x 10" 4 


Calculated Terminal Re Number (Eq. 13-51) 


Re v 




0.65-3.19 


0.15-0.85 


0.002-0.011 


Archimedes Number 


Ar 




83.89 


18.57 


0.192 


Particle Contact Duration 


tc 


s 


2.54 x 10~ 3 






Fluid CFL Condition 


tf 


s 


1.4-4.5 xl0~ 5 






Fluid-particle Relaxation Time 


td 


s 


7.47 x 10~ 2 


1.56 x 10" 3 


1.56 x 10" 4 



Also included in Tabled] are the relevant timescales for the simulations. The 
particle contact duration t c and fluid CFL condition t f are described in Sections 
13.21 and 13.11 respectively. The fluid-particle relaxation time is the characteristic 
time that a falling particle in Stokes flow will converge to its terminal velocity. 
This is given by td — m/b from Eq. (|34|) . This relaxation time provides another 
minimum timestep for the SPH-DEM simulation, given by 

AUelax < ^q"^- ( 36 ) 

The physical properties of the solid DEM particles are constant over all 
the simulated cases. Since the results of the test cases are insensitive to the 
particle-particle contacts, a relatively low spring stiffness of k = 10~ 4 kg/s 2 was 
used. This value ensures that the timestep is limited by the fluid CFL condition, 
rather than the DEM timestep, significantly speeding up the simulations. 

5. Single Particle Sedimentation (SPS) 

This section describes the results from SPH-DEM simulations using the first 
test case (SPS). We tested one and two-way coupling between the phases, the 
effect of different drag laws (Stokes and Di Felice), different fluid properties (air, 
water and water-glycerol) and the effect of varying the fluid resolution. 

5.1. One and two-way coupling in Stokes flow 

For a single particle falling in Stokes flow the standard Stokes drag equation, 
Eq. (|28p. can be used for the drag. Since Stokes drag law assumes a quiescent 
fluid, the force on the fluid due to the particle is set to zero (f = in Eq. ((26)) ). 
This implements a one-way coupling between the phases. Note that the SPH 
particles can still interact with the DEM particles through the porosity field, 
but for a single particle this effect will be negligible. 

In Figure 0] the evolution of a DEM particle's vertical speed in water is 
shown for one-way and two-way coupling. Also shown is the expected analytical 
prediction using Eq. (l34l) . The falling DEM particle reproduces the analytical 
velocity very well for both one-way and two-way coupling and the error between 
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Figure 4: Sedimentation velocity for a single particle in different fluids falling from rest with 
both one-way and two-way coupling. The dashed line is the theoretical result integrating 
Stokes law. The y-axis shows the particle vertical velocity scaled by the expected terminal 
velocity |ut| and the x-axis shows time scaled by the drag relaxation time t^. The inset shows 
the percentage error between the SPH-DEM and the expected trajectory. 



the two curves is less than 2% for the vast majority of the simulations. Note 
that the error curve does diverge to 5% when the particle is first released, but 
this is is a short-lived effect and the error drops below 2% after a time of about 
td-, the relaxation time for the drag force. 

These results indicate that the pressure gradient calculated from the SPH 
model, very accurately reproduces the buoyancy force on the particle, balancing 
out the drag force at the correct terminal velocity. The results are identical for 
both one-way and two-way coupling, indicating that the drag force on the fluid 
has a negligible effect. This is true as long as the fluid resolution is sufficiently 
larger than the DEM particle diameter, and this is explored in more detail in 
Section IQ 

Figure 2] also shows the same result for a DEM particle falling in air and 
in the water-glycerol mixture. For air, the drag force on the particle is much 
lower than for water, and the particles do not have time to reach their terminal 
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velocity before reaching the bottom boundary, where the simulation ends. As 
for the previous simulation with water, there is initially a larger (approx 4%) 
underestimation of the particle vertical speed, but once again this occurs only 
for a very small time period and does not affect the long term motion of the 
particle. For the majority of the simulation the error is less than 0.3% for both 
one-way and two-way coupling. 

The results for the water-glycerol fluid are qualitatively similar to water. 
Here the drag force on the particle is much higher than for water and the par- 
ticle reaches terminal velocity very quickly. However, as long as the simulation 
timestcp is modified to resolve the drag force relaxation time td as per Eq. (I36[) , 
the results are accurate. For both the one-way and two-way coupling, the sim- 
ulated velocity matches the analytical velocity very well and the error remains 
less than 0.3% for the duration of the simulation. 

In summary, the results for the one-way and two-way coupling between the 
fluid and particle for all the reference fluids are very accurate, and reproduce 
the analytical velocity curve within 0.3-2% error besides short-lived higher de- 
viations at the initial onset of motion. All data can scale using ut and td for 
velocity and time respectively. 

5.2. Effect of Fluid Resolution 

In this section we vary the fluid resolution to see its effects on the SPS results. 
Using water as the reference fluid, four different simulations were performed with 
the number of SPH particles was ranging from 10x10x15 particles to 40x40x60. 
Using the SPH smoothing length h as the resolution of the fluid, this gives a 
range of 1.5d < h < 6d, where d is the DEM particle diameter. 

Figure [5] shows the percentage difference between the average terminal ve- 
locity of the particle and the expected Stokes law. The error bars in this plot 
show one standard deviation of the fluctuations in the terminal velocity around 
the average. 

The h/d — 6 resolution corresponds to that used in the previous one and two- 
way coupled simulation, and the percentage error here is similar to the one-way 
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Figure 5: The effect of fluid resolution for the SPS test case, with water as the surrounding 
fluid. The x-axis is h/d, where h is the SPH resolution and d is the DEM particle diameter. 
The y-axis shows the average percentage error between the particle terminal velocity and the 
analytical value. The errorbars show one standard deviation from the mean. 

case, which is a mean of 0.2% with a standard deviation of 0.8%. As the fluid 
resolution is increased there is no clear trend in the average terminal velocity, 
but there is an obvious increase in the fluctuation of the terminal velocity around 
this mean. For h/d > 2, the standard deviation of this fluctuation is less than 
1%, but this quickly grows to 3% for h/d = 1.5. 

The increased error as the fluid resolution approaches the particle diameter 
is due to one of the main assumptions of the AVNS equations, i.e. that the fluid 
resolution length scale is sufficiently larger than the solid particle diameter. In 
this case the smoothing operator used to calculate the porosity field is also much 
greater than the particle diameter and this will result in a smooth porosity field. 
As the fluid resolution is reduced to the particle diameter the calculated porosity 
field will become less smooth and there will emerge local regions of high porosity 
at the locations of the DEM particles. As the fluctuations in the porosity field 
become greater this in turn will cause greater fluctuation in the forces on the 
SPH particles leading to a more noisy velocity field. 
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Another trend, that is not clear in these results but can be seen for solid 
particles with higher density, is the terminal velocity of the particle increasing 
with increasingly finer fluid resolution. Due to the two-way coupling, the drag 
force on the particle will be felt by the fluid as an equal and opposite force. This 
will accelerate the particles by a amount proportional to the relative mass of the 
SPH and DEM particles. For higher resolutions the mass of the SPH particles 
is lower, leading to an increase in vertical velocity of the affected fluid particles. 
Since the DEM particle's drag force depends on the velocity difference between 
the phases, which is now smaller, this will lead to a increase in the particle's 
terminal velocity. For the SPS test case shown here, the single particle does not 
exert too much force on the fluid and this is not a very large effect. As the fluid 
resolution is decreased from h/d = 6 to 2, there is a slight increase (on the order 
of 1-2%) in the terminal velocity, but lower than this the trend is lost, likely due 
to the increasingly noisy results due to the fluctuations in the porosity field. 

5. 3. The effect of fluid properties and particle Reynolds number 

We have used three different reference fluids in the simulations, correspond- 
ing to air, water and a water-glycerol mixture. Using the SPS test case, this 
results in a range of particle Reynolds numbers between 0.011 (water-glycerol) 
and 3.19 (air), allowing us to explore a realistic range of particle Reynolds 
numbers. We have further extended this range by considering two additional 
(artificial) fluids with a density of water but lower viscosities, resulting in a 
range of 0.011 < Re p < 9. 

Rather than assuming Stokes flow as in the previous sections, here we 
will use the Di Felice drag law (e = 1), which is assumed to be valid for all 
Reynolds numbers. This will be compared against fully resolved simulations 
using COMSOL Multiphysics (finite element analysis, solver and simulation 
software, http : //www . comsol . com/) . 

Figure [6] shows the average error in the terminal velocity measured from the 
SPH-DEM simulations using both the Stokes and Di Felice drag laws, using the 
COMSOL results as the reference terminal velocity. Since the two drag laws 
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Figure 6: Error in SPH-DEM average SPS terminal velocity at different terminal Re numbers. 
The fully resolved COMSOL simulation is used as reference for the error calculation. The 
solid red and dashed green lines show the results using either Stokes or Di Felice drag law. 
The dotted blue line shows the reference terminal ve locity calculated using the Coulson and 
Richardson drag law (Coulson and Richard son! . 1 199*3) . 

are equivalent at low Re p , they give the same result at Re p = 0.01. As Re p 
increases, the plots diverge, and the simulated terminal velocity using the Stokes 
drag quickly becomes much larger than the COMSOL prediction (as expected 
since the Stokes drag law is only valid for low Re p ). In contrast, the Di Felice 
drag law results in a simulated terminal velocity that follows the same trend as 
the COMSOL results. At low Re p the DEM particle falls slightly (~ 5%) faster, 
at higher Re p it falls slightly (3-6%) slower. 

For further comparison, the COMSOL results have also been compared with 
the analytical drag force model proposed by Coulson and Richardson and re- 
produced in Eq. (|29j) . The expected terminal velocity was calculated using this 
model and plotted alongside the SPH-DEM results in Figure [6] As shown, the 
COMSOL results agree with this analytical terminal velocity to within 3.5% 
over the range of Re p considered. 

While the results in previous SPS sections have shown the SPH-DEM model 
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can accurately reproduce the expected terminal velocity assuming a given drag 
law (Stokes), this section illustrates that the final accuracy is still largely de- 
termined by the suitability of the underlying drag law chosen. However, a full 
comparison of the numerous drag laws currently in the literature is beyond the 
scope of this paper, and for the purposes of validating the SPH-DEM model we 
can assume that the chosen drag law (from here on the Di Felice), approximates 
well the true drag on the particles. 

6. Sedimentation of a Constant Porosity Block (CPB) 

This section shows the results from the Constant Porosity Block (CPB) test 
case. In a similar fashion to the SPS case, we explore the effect of fluid resolution 
and fluid properties. In addition, we consider the influence of a new parameter, 
the porosity of the block, on the results. All the simulations in this section use 
two-way coupling, as the hindered fluid flow due to the presence of the solid 
particles is an important component of the simulation. As the porous block 
falls, the fluid will be displaced and flow upward through the block, affecting 
the terminal velocity. All the simulations use the Di Felice drag law, which is 
necessary to incorporate the effects of neighbouring particles (lower porosity) 
on the drag force. 

Figure [7] shows an example visualisation during the simulation of a block 
with porosity e = 0.8 falling in water. On the left hand side of the image are 
the DEM particles (coloured by porosity e») falling in the fluid column. The 
porosity of most of the DEM particles is e = 0.8, as expected, except near the 
edge of the block where the discontinuity in particle distribution is smoothed 
out by the kernel (with smoothing length h c = 6d) in Eq. ([9]). This results in 
a porosity greater than 0.8 for DEM particles whose distance is lower than h c 
from the edge of the block. We will show in Section 16.11 that this effect can be 
limited/avoided by choosing a smaller smoothing length. 

On the right hand side is shown a vector plot of the velocity field at x — 0. 
This shows the upward flow of fluid due to the displacement of fluid by the 
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Figure 7: Visualisation of the DEM particles for the Constant Porosity Block test case. On the 
left the DEM particles are shown coloured by porosity £j, and a transparent box representing 
the simulation domain. On the right the corresponding fluid velocity field is shown at x = 0, 
with the arrows scaled and coloured by velocity magnitude. 
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Figure 8: Average percentage change in the terminal velocity and average porosity of the 
Constant Porosity Block (CPB) with e = 0.8 in water for varying fluid resolution. Errorbars 
in the terminal velocity points show one standard deviation of the vertical velocity data from 
the average, taken over a time period of 0.34 s after the terminal velocity has been reached. 

particles as they fall. Also noticeable are fluctuations in velocity near the edges 
of the block, which are discussed in more detail in Section [6.31 

Shortly after release, the vertical velocity of the CPB converges to a termi- 
nal velocity that is consistent with the expected terminal velocity, although it 
is slightly (less than 5%) higher than expected. The systematically increased 
terminal velocity is due to reduced drag at the edges of the block due to the 
finite width of the smoothing kernel. As the width of the smoothing kernel h 
used to calculate the porosity field is larger (by a factor of 2-6, see Figure [8] for 
details) than the particle diameter d, the porosity field near the edges of the 
CPB will be smoothed out according to the width of the kernel. This results in 
a slightly higher apparent local porosity and a reduced drag than what would 
be expected with e = 0.8. 

6.1. The effect of fluid resolution 

Figure [8] shows the percentage difference between the vertical velocity of 
the block and the expected terminal velocity. The results from five different 
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simulations are shown, each with a different fluid resolution ranging from h/d = 
6 to h/d = 2. The porosity is set to e = 0.8. The h/d — 6 plot uses a 
smoothing kernel that is 6 times greater than the DEM particle diameter, leading 
to increasing smoothing of porosity field near the edges of the block. Integrated 
the porosity field over the volume of the CPB leads to a porosity of 0.85, about 
6% higher than the true porosity of the block. This results in an increase of 
22% in the terminal velocity of the block. Increasing the fluid resolution to 
h/d = 5 causes the error to decrease to 15%, since the interpolated porosity at 
the edge of the block is now closer to the set value of e = 0.8. Further increases 
in the fluid resolution consistently decreases the measured terminal velocity 
until at h/d = 2 the error is only 5% of the expected value. These results 
illustrate how the smoothing applied to the porosity field can have dramatic 
results on the accuracy of the simulations. This is largely due to the fact that 
the modelled drag only depends on the local porosity, which does not properly 
consider the influence of porosity gradients on the applied drag force. Therefore 
the accuracy of the drag law near large changes in porosity is highly dependent 
on the magnitude of smoothing applied to the porosity field. This is true for the 
Di Felice law and the most oth er drag laws pro posed in the literature, but there 



has been some recent work by 



Xu et al 



(2007), which attempts to account for 



the influence of the porosity gradient, but we will not study this further here. 
6.2. The effect of porosity 

Varying the porosity of the CPB allows us to evaluate the accuracy of the 
SPH-DEM model at different porosities. Figure [5] shows the average terminal 
velocity of the block, as measured from SPH-DEM simulation of the CPB over 
a range of porosities from e = 0.6 to 1.0. Results using both water and water- 
glycerol as the interstitial fluid are shown on the same plot by scaling the y-axis 
by the expected terminal velocity of a single DEM particle. The average terminal 
velocity is taken after the block has reached a steady terminal velocity and the 
error bars show one standard deviation of the vertical velocity from the average. 
For these simulations, again h/d = 2 is used. 
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Figure 9: Average terminal velocity (scaled by |ut[, the expected terminal velocity of a single 
DEM particle) of the Constant Porosity Block (CPB) in water and water-glycerol for varying 
porosity and h/d = 2. Error bars show one standard deviation of the vertical velocity data 
from the average, taken over a time period of 0.34 s. The y-axis is scaled by |ut|, the expected 
terminal velocity of a single DEM particle given by Eq. (1351) . which corresponds to the SPS 
test case. 
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Shown with the SPH-DEM results is a reference line showing the expected 
terminal velocity computed using Eq. (|35[) and the input porosity of the block. 
The SPH-DEM results for both water and water-glycerol match this reference 
line very well over the range of porosities tested. At lower porosities the vertical 
velocity of the CPB suffers from increasing fluctuation around the mean. This 
is a consequence of fluctuations seen in the surrounding fluid velocity, and will 
be described further in Section IQ1 

In summary, the simulated terminal velocity for the CPB matched the ex- 
pected value over the range of resolutions and porosities considered, as long as 
the resolution of the fluid phase (give by h) is sufficient to resolve the porosity 
field of the given problem. For the CPB we have a jump at the edge of the block 
from the given porosity of the block to the surrounding e = 1. We found that 
as long as the fluid resolution was kept at h = 2d, where d is the DEM particle 
diameter (i.e.. the length scale of the porosity jump), the results matched the 
theoretical predictions within 5% error. 

6.3. Effect of Porosity Gradients on Fluid Solution 

In the previous section it was shown how the smoothing of the porosity 
discontinuity of the block slightly affected the drag on the DEM particles and 
the final terminal velocity of the block. In this section we will show how the 
high porosity gradients near the edge of the block also gave rise to further effects 
on the SPH solution for the fluid. 

Figure [10] shows the vertical velocity and porosity for all the SPH particles in 
a CPB simulation with fluid resolution h/d = 2 and porosity e = 0.8. These val- 
ues are plotted against the vertical position of the SPH particles. The porosity 
is rather smooth and clearly shows the location of the CPB. However, there are 
fluctuations in the vertical velocity of the SPH particles near the edges of the 
block, much larger than the rather small average (positive) velocity inside the 
block. These fluctuations are present to different degrees in all of the SPH-DEM 
simulations and their magnitude is proportional to the local porosity gradient. 
Therefore, their effect is strongest for the simulations with low porosity or fine 
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Figure 10: Scatter-plot of the vertical velocity (red dots) and porosity (green line) versus 
height for all the SPH particles. The test case was CPB with a porosity of e = 0.8 in water 
as the surrounding fluid, the fluid resolution was h/d = 2 and a ar t =0.1 



fluid resolution (i.e. small h c ). 

Given the correlation of these fluctuations with high porosity gradients, their 
sourc e is likely to be due to errors in the SPH pressure field. It is well-known 
(e.g. (|Colagrossi and Landrinil . 120031 )) that SPH solutions can exhibit spurious 
fluctuations in the pressure field, which normally have little or no effect on 
the fluid velocity. For our simulations the pressure of each SPH particle is 
proportional to (p/epo) 7 and is therefore very sensitive to changes in e. It is 
likely that for high porosity gradients the pressure variations that are normally 
present would be amplified and generate corresponding large fluctuations in the 
velocity field. 

As long as the fluctuations do not grow too large, they do not effect the mean 
flow of the fluid, as evidenced by the accurate reproduction of the expected 
terminal velocity in the previous sections. To ensure the simulation accuracy, it 
was found that the application of an artificial viscosity with strength a ar t =0.1, 
see Eq. (fl"3"|) . was enough to damp out the fluctuations in velocity so that they 
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did not have a significant effect on the results. This value of a ar t was used in all 
of the CPB simulations shown here. The artificial viscosity has no effect on the 
settling velocity of the SPS or CPB since this viscosity is only applied between 
SPH particles and is not included in the fluid-particle coupling term (Eq. I!?3"|) . 
However, for systems where the fluid viscosity plays an important role (e.g. the 
Rayleigh Taylor instability in the following section), this has an effect which 
will be described in the next section. 

7. Rayleigh- Taylor Instability (RTI) 

The classic Rayleigh- Taylor fluid instability is seen when a dense fluid is 
accelerated into a less dense fluid, for example, under the action of gravity. 
Consider a water column of height h filled with a dense fluid with density pd and 
viscosity Vd located above a lighter fluid with parameters pf and Vf. For the RTI 
test case, the lower and higher density fluids are represented by the pure fluid 
and the suspension, respectively. If the height of the interface between the two 
fluids is perturbed by a normal mode disturbance with a certain wave number 
k (see Figure QT] and Eq. (J39J) ) , then this disturbance will grow exponentially 
with time. 

The two-fluid model of a Rayle i gh- Ta ylor instability was derived in the au- 



thoritative text by IChandrasekharl (|196ll ). The exponential growth rate n(k) of 



a normal mode disturbance with wave number k at the interface between the 
two fluids (with zero s urface tension) is characterised by the dispersion relation 



rtChandrasekhar 



19611 ) given by 



^( a f - a d) + 1 



{0L c q d + a f q c — k) - 4ka f a d 



4fc 2 

+ ——{a f v f - u d v d )[a d qf - a f q d + k(a f - a d )} 
4fc 3 

+ -^{ a S v i - a d Vd) 2 {qf ~ k)(q d - k) = 0, (37) 



where v^d — fJ>f,d/ Pf,d is the kinematic viscosity of the two phases, a/. 
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Figure 11: Diagram showing a cross-section of the initial setup for the Rayleigh- Taylor In- 
stability (RTI) test case. The upper grey area is the particle-fluid suspension with effective 
density and viscosity and fi^ , the lower white region is clear fluid with density and viscosity 
p j and /if . The suspension is given an initial vertical perturbation with wave number k and 
amplitude d/i. 



Pf,d/(Pf + Pd) is a density ratio and qj d = k 2 + n/isfd is a convenient abbrevi- 
ation. 

For this test case, we use an identical initial condition as in the CPB test 
case, with an initial block of particles immersed in the fluid with an initial 
porosity of e = 0.8. Using the density of the surrounding fluid p/, the effective 
density of the fluid-particle suspension is pd = epf + (1 — e )Pp- This system can 
be approximated using a two fluid model, where the suspension is treated as a 
fluid with density pd and kinematic viscosity v&. Therefore, in a similar fashion 
to the RT instability described above, it is expected that an initial disturbance of 
the interface between the two "fluids" will increase with an exponential growth 
rate which is a solution of Eq. (|37|) . 

The effective visc osity of the sus pension p.^ is estimated here using Krieger's 
hard sphere model (|Kriegeri 119591 ) (assumed to be valid for both dilute and 
dense suspensions) 
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Figure 12: Visualisation of the DEM particles (left) and velocity field (right) at x = in the 
x-z plane, for the Rayleigh Taylor (RT) test case, using e = 0.8 and water-glycerol as the 
surrounding fluid. The growth rate for this simulations versus time can be seen in Figure [14] 



, 2.5(l-e m ;„) 



= M/ — , (38) 

where e m j„ = 0.37 is the porosity at the maximum packing of the solid 
particles. 

We generate an initial disturbance in the interface between the two "fluids" 
by adding a small perturbation to the vertical position of every DEM particle 

Azi = -^(1 - cos(fe K a; i ))(l - cos(k y yi)), (39) 

where k x — k y = 2tt/w and Xi and yi are the x- and y-coordinates of the 
position of particle i. This yields a symmetric disturbance in the interface with 
a wave length equal to the box width w. 
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Figure 13: Growth of Rayleigh- Taylor instability using water. The red pluses and green 
crosses show the position of the lowest DEM particle when the artificial viscosity is either 
added or not. The two reference lines show the expected growth rate using the lowest and 
highest porosity of the CPB. 



Figure [T2l shows the positions of the DEM particles during the growth of the 
instability, along with the fluid velocity field at x = 0. At this time there is a 
strong fluid circulation that is moving downward in the centre of the domain 
and upward at the edges (not visible in this cut). This causes the growth of 
the instability by increasing the sedimentation speed of the DEM particles near 
the centre while reducing or even reversing the sedimentation of those particles 
near the outer boundaries of the domain. The movement of the DEM particles 
matches the expected behaviour of the instability, and the wave length of the 
dominant mode is identical to the initial perturbation given to the DEM particle 
positions. Next we will attempt to quantitatively compare the SPH-DEM results 
to the growth rate predicted by the analytical two-fluid model. 

In Figure IT3l the growth of the RT instability versus time for e = 0.8, fluid 
resolution h/d = 2 is shown using water as the surrounding fluid. The symbols 
show the vertical position of the lowest DEM particle, which provides an ap- 
proximate measure of the instability amplitude. The vertical displacement of 
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this point over time can be compared with the estimated growth rate for the 
RT instability as given by the two- fluid model in Eq. (|3"T1) . The growth rate 
of the instability is added to the expected sedimentation speed using Eq. (|35p 
to calculate the expected trajectory of the lowest DEM particle. Using the pa- 
rameters of the simulation (including e = 0.8) and solving for the growth rate 
leads to a growth curve given by the lowest blue dashed line. The two-fluid 
model is included here as a benchmark, but it should be noted that this model 
contains some significant approximations in treating the particle suspension as 
an equivalent fluid, and is not necessarily more accurate than the SPH-DEM 
results. While a constant porosity of 0.8 is used for the two-fluid RTI model, 
the porosity of the DEM particles ranges from 0.8 < e < 0.86 at t = (ini- 
tial conditions) and the porosity at the leading front of the instability grows 
over time, reaching a value of 0.93 at the timestep shown in Figure [T^] and a 
maximum value of 0.95 before the instability meets the bottom boundary. To 
account for this variation in porosity, we instead use the analytical model to 
obtain an upper and lower bound to the instability growth. The upper bound 
is calculated using epsilon — 0.8 (the blue dashed line) and the upper bound is 
calculated using e = 0.93, which gives the purple dashed line. 

The SPH-DEM results are shown for the cases where the artificial viscosity 
is either applied (a art = 0.1) or not used (a art = 0.0). In both cases there is a 
clear exponential growth of the RT instability and only the quantitative growth 
rate differs between the two simulations. If the artificial viscosity is applied, the 
growth rate of the instability is lower than both of the two reference bounds. If 
the artificial viscosity is not used, the growth rate of the instability is increased 
to lie between the two bounds. After t = 0.15 s the growth rate becomes slower 
than the upper bound, but by this time the bottom of the instability is close to 
the bottom boundary, and we do not expect the two-fluid model (which assumes 
an unbounded domain) to apply. 

Figure [TJ] shows the same results but using water-glycerol as the interstitial 
fluid. In this case the physical viscosity of the fluid is proportionally greater 
than the artificial viscosity applied, and therefore the addition of the artificial 
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Figure 14: Growth of Rayleigh- Taylor instability using water-glycerol. The red pluses and 
green crosses show the position of the lowest DEM particle when the artificial viscosity is 
either added or not. The two reference lines show the expected growth rate using the lowest 
and highest porosity of the CPB. 



viscosity has a lesser effect. For both a ar t = 0.1 and a ar t = 0.0 the growth rate 
of the instability lies between the two bounds, except when the DEM particles 
reach the bottom of the domain and wall effects start to dominate. 

While it is encouraging that the SPH-DEM results (without artificial viscos- 
ity) closely match the expected growth of the RT instability, the results highlight 
the negative effect of the artificial viscosity when used in problems where the 
fluid viscosity is important. Therefore it is necessary to develop an alternative 
method of reducing the velocity fluctuations near high porosity gradients, and 
this is the subject of current work. However, it is important to note that for 
the majority of applications the addition of a small amount of artificial viscos- 
ity has no significant effect on the results and is successful in eliminating the 
problematic velocity fluctuations. 

In summary, the results from the RTI simulations using water-glycerol show 
that the SPH-DEM simulation can accurately reproduce the Rayleigh- Taylor 
instability. The addition of an artificial viscosity, while successful in dampen- 
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ing the velocity fluctuations, increases the effective viscosity of the system and 
reduces the growth rate of the instability. 



8. Conclusion 

We have presented an SPH implementation of the locally averaged Navier 
Stokes equations and coupled this with a DEM model in order to provide a 
simulation tool for one or two-way coupled fluid-particle systems. One notable 
property of the resulting method is that it avoids the use of a mesh and is 
completely particle-based. It is therefore suitable for those applications where a 
mesh presents additional problems, for example, fr ee surface flow or flow around 
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complex, moving and/or intermeshed geometries ([Robinson et al. 

The SPH-DEM formulation was applied to 3D single and multiple particle 
sedimentation problems and compared against analytical solutions. For single 
particle sedimentation the SPH-DEM simulations reproduced the analytical so- 
lutions very well, with less than 2% error over a wide range of Particle Reynolds 
Number 0.011 < Re p < 9 and fluid resolutions. Only when the fluid resolution 
became less than or equal to 1.5 times the particle diameter did the results start 
to diverge from the expected solution. 

For the multiple particle sedimentation test case using the Constant Porosity 
Block (CPB), the SPH-DEM method accurately reproduced the expected ter- 
minal velocity of the block, over a range of porosities 0.5 < e < 1.0 and Particle 
Reynolds Number 0.002 < Re p < 0.85. This is a general consequence of the 
locally averaged Navier Stokes equation, and is not specific to the SPH-DEM 
method. The minimum resolution of the porosity field will always be much 
coarser than the DEM particle diameter, so any discontinuities in the particle 
distribution will be smoothed and thus not be accurately reflected in the poros- 
ity field. Furthermore, the results from the CPB also showed an instability in 
the SPH fluid phase that occurred near the edges of the block. The high porosity 
gradients in this region give rise to fluctuations in velocity of the SPH particles, 
which are likely due to fluctuations in the pressure field being amplified by the 
sudden change in porosity. Adding a small amount of artificial viscosity to the 
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simulations was sufficient to damp these fluctuations and prevent them from 
affecting the terminal velocity of the block. 

The Rayleigh- Taylor Instability (RTI) test case successfully reproduced the 
instability and its growth rate for both water and water-glycerol. For this test 
case the addition of artificial viscosity was not necessary, due to the relatively 
high porosity e = 0.8 and lower porosity gradients at the interface between the 
suspension and clear fluid. 

Removing the SPH velocity fluctuations near high porosity gradients is the 
subject of current work, and promising results have already been obtained by 
either calculating the drag separately on the fluid or re-deriving the SPH equa- 
tions from a Lagrangian formulation. Besides this issue, it was found that the 
SPH-DEM model successfully reproduced most of the expected results from the 
analytical test cases over a wide range of Reynolds Numbers and porosities, and 
promises to be a flexible and accurate tool for modelling particle-fluid systems. 

In the future, the me thod will be applied t o dispersion of solids in fluid 



or fluid-gas environments ([Robinson et al 



20121 ). Other relevant directions for 



future developments are: choice of drag law and the inclusion of the added mass 
and lift forces; the choice of DEM particle contact forces and the inclusion of 
friction and lubrication forces; and the inclusion of surface tension effects. These 
questions require further study of the method and the choice of the parameters, 
laws and assumptions. 
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